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In chaotic dynamical systems, a number of rare trajectories with low level of chaoticity are embedded in chaotic sea, 
while extraordinary unstable trajectories can exist even in weakly chaotic regions. In this study, a quantitative method 
for searching these rare trajectories is developed; the method is based on multicanonical Monte Carlo and can estimate 
the probability of initial conditions that lead to trajectory fragments of a given level of chaoticity. The proposed method 
is successfully tested with four-dimensional coupled standard maps, where probabilities as small as 10 -14 are estimated. 
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1. Introduction 

In chaotic dynamical systems, a number of rare trajec- 
tories with low level of chaoticity are embedded in chaotic 
sea, which show apparently small values of Lyapunov ex- 
ponents. Also, there can be extraordinary unstable trajec- 
tories with apparently large values of Lyapunov exponent 
even in weakly chaotic systems. 

Numerical search for such rare trajectories and quan- 
titative estimation of their probabilities are important for 
the understanding of these dynamical systems. Even when 
Lyapunov exponent is converged to a unique value almost 
everywhere in a chaotic sea, the probability of finite- length 
trajectories of a given chaoticity reflects transient behav- 
ior of the system and hence provides information on fine 
structures in the state space [13|. 

In this paper, we develop a method for calculating 
probabilities of initial conditions that lead to trajectory 
fragments of a given level of chaoticity. This method is 
based on the multicanonical Monte Carlo algorithm 
which is a version of dynamic Monte Carlo (Markov 
chain Monte Carlo) [H, [jj and provides a powerful tool 
for calculating small probabilities under a given probabil- 
ity measure. The proposed method is tested with four- 
dimensional coupled standard maps. It is shown that the 
method can deal with rare events of very small probabili- 
ties such as ~ 10 -14 . 

The method proposed in this paper can be regarded as 
a variant of the method recently introduced in Yanagita 
and Iba [5j; the proposed method uses the multicanon- 
ical Monte Carlo instead of the replica exchange Monte 
Carlo used in [5J. A novelty in this paper is that proba- 
bilities of initial conditions that generate rare trajectories 



are successfully estimated, while no quantitative result is 
obtained in the previous paper [Hj]. Also, the multicanon- 
ical algorithm has advantages over the replica exchange 
Monte Carlo; it enables direct computation of the desired 
probabilities and is efficient in cases with first order tran- 
sitions [l[. 

Our studies are partly motivated by the studies Q, 
where rare structures with high or low chaoticity are ex- 
plored by simulating fictitious particles that are moved, 
split, and systematically erased; their method can be re- 
garded as tracking rare "pseudo-orbits" in a parallel man- 
ner. The method proposed in this paper samples initial 
conditions and seems complementary to their approach. 
There are also references [JS Ell [III that discuss sampling 
of unstable periodic orbits in chaotic systems by dynamic 
Monte Carlo or related methods. None of these studies, 
however, seems to compute probabilities of rare trajectory 
fragments using the multicanonical algorithm. 

2. Algorithm 

2.1. Forgetting Time as a Measure of Chaoticity 

Let us consider a deterministic dynamical system 
x(t + 1) = ip(x(t)) with discrete time t = 1,2,... and as- 
sume a trajectory x(t) originated from an initial condition 
x(0) — xq. Specify a measure / for the chaoticity of the 
trajectory x{t) and consider it as an integer- valued func- 
tion /(xq) of the initial condition Xq\ for a real-valued /, 
appropriate discretization of its value is assumed. Then, 
our interest is in estimating the probability 
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where 6 is defined as S(s) = 1 if s = 0; S(s) — otherwise. 
dxo is uniform measure on the space of initial conditions 
and D is the volume of the entire space. 
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An important point is the choice of /(xo); assuming fi- 
nite precision arithmetic of computers, it should be some- 
thing like "approximate Lyapunov exponent" of a finite 
piece of the trajectory starting from xq- 

A possible way is to fix the length T of the trajecto- 
ries and define f(xo) as the largest eigenvalue At(xo) of 
the product Ylt=i °f tne Jacobian matrices 3(x(t)) 

along the piece of the trajectory x(t) starting from xq. 
Here (i, j)-component of J(x(i)) is given by Jy(x(i)) = 

^ , where tpi and ith component of the 

0x i x=x(i) 

function ip and its argument, respectively. This way of 
choosing f(xo) is, however, not suited for numerical com- 
putation when Xt(xq) is strongly dependent on xq. When 
Xt(xo) is large we should choose small T to control the 
effect of round-off error. On the other hand, larger T is 
desirable in the region of small At(xo) to filter out trajec- 
tories with low chaoticity. Thus, it is difficult to choose 
the value of T adequate for all values of At(xo). 

Our solution is to reverse the idea and use the min- 
imum value T e (xo) of T that satisfies Xt(xo)c > 1 as a 
measure of chaoticity, where e is a small constant beyond 
the machine epsilon; larger value of T e (xo) means larger 
number of time-steps required to forget the initial con- 
dition and implies lower chaoticity. The switching from 
f(x ) = X T (x ) to "forgetting time" f(x ) = T e (x ) pro- 
vides a computationally stable criterion of chaoticity, be- 
cause by definition the latter is insensitive to round-off 
error. 

2.2. Multicanonical Algorithm 

We introduce the multicanonical algorithm to estimate 
the probability P(/); the algorithm consists of two stages, 
training and measurement. In the training phase, we 
construct the approximate probability P(f) step-by-step 
through dynamic Monte Carlo simulations. At each step, 
we perform Monte Carlo simulation with the weight func- 
tion l/P(/(xo)) using the current estimate of P(/). If 
P{f) is a good approximation to P(/) in a prescribed in- 
terval of /, the histogram h(f) is almost flat in the interval, 
because 

Kf) i / S(f(x ) - f)^z dx = Oft ~ 1. 

If h(f) is not sufficiently flat, we modify P(/) until an 
almost flat histogram of / is obtained. Several methods are 
proposed for efficient tuning of the weight in the training 
phase; in this study, we use a method due to Wang and 
Landau [l2T |. After we obtain a good approximation of 
P(/) we enter the measurement phase and perform a long 
run of simulation with a fixed weight 1/P(/(xq)). Then, 
the final estimate P*(f) of P(f) is obtained as P*(f) oc 

Kf)P(f). 

The key to the implementation of the algorithm is that 
the weight l/P(/(xo)) of xq is expressed with a compos- 
ite function of a univariate function P(f) and a known 



function /(xo); it enables easy adaptation of the weight 
using outputs of the simulations. On the other hand, the 
multicanonical algorithm enjoys fast mixing of the Markov 
chain; sampling in a wide range of / realizes a kind of "an- 
nealing" effect and the convergence becomes much better 
compared with the cases where only the tail regions are 
sampled. 

So far we discuss the generic algorithm. To implement 
the multicanonical algorithm in the present case, where 
/(xo) = T e (xo), we should specify dynamic Monte Carlo 
algorithm for the sampling with the weight l/P(T e (xo)). 
Here we use the Metropolis algorithm @, in which a 
candidate Xq cw is generated using a hierarchical proposal 
distribution used in The candidate Xq £W is accepted 
if and only if a uniform random number r € [0, 1) satis- 
fies r < £f~'f X ° c J), , where xg w is the current value of x ; 

to compute T e (xQ ew ) at each step of the Metropolis algo- 
rithm, we simulate the trajectory from the initial condition 
x% ew until X T {x^ ew )e > 1. 

2. 3. Computation of the Largest Eigenvalue 

To implement the proposed algorithm for high- 
dimensional dynamical systems, we should compute the 
largest eigenvalue A of the matrix 3t — Y[J=i em ~ 
ciently. In this study, we approximate it using the power 
method as A ~ | J™£| |, where £ is a constant vector of unit 
length ||£|| = 1. We found that m = 1 often gives a good 
approximation. 

3. Numerical Experiments 

3.1. Searching Low Chaoticity 

Let us consider four-dimensional coupled standard 
maps 

K b , . . 

u n+ x = u n - — sm(27ru„) + — sm(27r(w„ + y n )), 
Zir Ztt 

Vn+l = V n + U n+ i, 

K . ,„ s b . , n . 
x n +i = x n - — sm(27ry„) + — sm(27r(u„ + y n )), 

Z7T Z7T 

Un+1 = Vn + X„+l, 

which is a well-studied example of volume preserving 
maps. 

First, we test the efficiency of the proposed algorithm 
to find tiny tori in chaotic sea. In Fig. [TJ a pair of tori in 
chaos found by the algorithm is plotted on the (u n ,v„)- 
plane; values of parameters K — 7.8 and b = 0.001 are 
chosen that most of the state space is covered by a chaotic 
sea. The threshold e is 2~ 43 . The values of T e correspond- 
ing to the initial conditions that lead to these tori are 
larger than 199 and the probabilities P(T £ ) are as small 
as 10~ 12 . The number of initial conditions tested in the 
proposed method is about 4 x 10 9 . 
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Figure 1: A pair of tiny tori in a chaotic sea found by the proposed 
method. Projections on the (u n , v n )-plane are shown. Enlargement 
of a tiny area in the small circle in the left panel is given in the right 
panel; further enlargement is given in the lower panel. K = 7.8 and 
b = 0.001. 



3.2. Probabilities P(T e ) 

An advantage of the multicanonical algorithm is that 
it computes probabilities P{T e ) of T e under uniform sam- 
pling of initial conditions. The probability of trajectory 
fragments that stays long time near tiny tori and/or can- 
tori will reflect fine structures of phase space, which can 
be quantified by P(T e ). In Fig.H P(T e ) for the model (JTJ 
is plotted. The result indicates that, in this parameter 
range, the volume of small regular regions decrease under 
increasing K, while increase under increasing b. 

The number of initial conditions tested in the proposed 
algorithm is 3 x 10 9 ~ 3 x 10 10 for each curve including 
training phase and rejected candidates; the corresponding 
computational times are 17 ~ 125 hours using a single core 
of Intel Xeon X5365, which is a quad-core CPU. These re- 
sults show that the proposed method can compute proba- 
bilities down to ~ 10~ 14 within reasonable computational 
times. 

For the comparisons, results of naive random sampling 
are also shown in Fig. [5] For (K,b) = (6.0,0.1), there is 
an overall agreement of the outputs following from the two 
sampling approaches. For (K,b) = (7.8,0.1), both results 
are also consistent, but naive random sampling with com- 
parable computational efforts (1.4 x 10 10 initial configura- 
tions) gives meaningful results only in a high probability 
region, as seen in the upper panel of Fig. [5] 

In these examples, we set the number m of iteration 
defined in Sec. 12.31 to unity; m = 2 and 3 are also tested 
and slight differences are found in some cases. 

4. Summary and Discussion 

A quantitative method based on multicanonical Monte 
Carlo is proposed for searching rare trajectories in chaos. 
The proposed method is tested with four-dimensional cou- 
pled standard maps and successfully computes the proba- 
bility of the forgetting time T e down to ~ 10~ 14 . 



Applications of the proposed method to dissipative 
and/or multi-basin systems will be interesting, as well as 
search for highly unstable trajectories in weakly chaotic 
systems. Another interesting subject is to develop a way 
to relate P(T e ) to "escape rate" and hence interpret it 
in the thermodynamic formalism [l3j ]. which connects our 
approach to existing studies on large deviations in chaos. 
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Figure 2: The probability P(T e ) is plotted as a function of A"(upper 
panel, b = 0.1) and fe(lower panel, K = 7.8), where e = 2~ 43 . The 
results of naive random sampling are shown by the symbol +, only 
when more than 10 samples are available. The symbols • at the right 
edge of the plot indicate values of the probability P(T f > 200). 
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